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Abstract 

The flavor composition of ultra-high energy cosmic neutrinos (UHECN) carries precious informa- 
tion about the physical properties of their sources, the nature of neutrino oscillations and possible 
exotic physics involved during the propagation. Since UHECN with different incoming directions 
would propagate through different amounts of matter in Earth and since different flavors of charged 
leptons produced in the neutrino-nucleon charged- current (CC) interaction would have different 
energy-loss behaviors in the medium, measurement of the angular distribution of incoming events 
by a neutrino observatory can in principle be employed to help determine the UHECN flavor ratio. 
In this paper we report on our investigation of the feasibility of such an attempt. Simulations were 
performed, where the detector configuration was based on the proposed Askaryan Radio Array 
(ARA) Observatory at the South Pole, to investigate the expected event-direction distribution for 
each flavor. Assuming f^-z^r symmetry and invoking the standard oscillation and the neutrino 
decay scenarios, the probability distribution functions (PDF) of the event directions are utilized to 
extract the flavor ratio of cosmogenic neutrinos on Earth. The simulation results are summarized 
in terms of the probability of flavor ratio extraction and resolution as functions of the number of 
observed events and the angular resolution of neutrino directions. We show that it is feasible to 
constrain the UHECN flavor ratio using the proposed ARA Observatory. 
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I. INTRODUCTION 



The observed energy spectrum of comic rays has been extended to beyond lO^^eV P[2], 
but httle is known about their origins and acceleration mechanism, which are important 
questions in astrophysics [3J. Ultra-high energy cosmic rays (UHECRs) are thought to be of 
extragalactic origin, such as being produced by active galactic nuclei (AGNs) or gamma-ray 
bursts (GRBs) [4J. Such UHECRs can generate ultra-high energy {E > 10^ ''eV) neutrinos 
via photo-pion production or proton-proton interaction: 

p-F7-)-A+-)-n-F 7r+, 
p + p ^ n^n^ir^, 

and the subsequent decays of charged pion and muon, e.g., 

where the targets can be the intergalactic medium near the astrophysical sources [51 E] (see 
Ref. [7j for a review), or the cosmic microwave background (CMB) photons (the Greisen- 
Zatsepin-Kuzmin process) [U [9]. Neutrinos originating from the GZK process are known 
as the cosmogenic neutrinos (or GZK neutrinos), which are guaranteed to exist based on 
the fact that both initial-state particles, i.e., the UHECR and the CMB photon, have been 
observed and that the notion is consistent with the observed GZK cutoff in the cosmic ray 
spectrum pTO]. Since the production of ultra high energy cosmic neutrinos (UHECNs) are 
tightly connected with UHECRs, such neutrino spectrum can help to resolve the puzzles of 
cosmic rays such as their composition pTl [T2] , the energy spectrum at the sources, and the 
cosmological evolution of the sources [131 IE] • 

Besides the overall spectrum, the relative flux ratio between different neutrino flavors, or 
briefly, the flavor ratio, can also provide information about the physical properties of UHECR 
sources. For example, the transition of flavor ratio at the source from : : = 1 : 2 : 
(pion source) to : 1 : (muon-damped source) with increasing energies due to synchrotron 
energy loss of muons can be used to constrain the strength of cosmic magnetic field [TSHTT] 
(neutrinos and antineutrinos are counted together because they are hard to be discriminated 
in the UHE neutrino detection). Furthermore, during the propagation from the source to the 
Earth, the flavor composition of a neutrino would oscillate [18] , and may even be altered by 
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some new physics beyond the Standard Model, such as the neutrino decay [T91423] . pseudo- 
Dirac states of neutrinos [21], sterile neutrinos with tiny mass differences [23], the violation 
of CPT or Lorentz invariance [201 ESI [2S] , and the quantum decoherence [231 [2Sl [2Z] (see 
Ref. [28] for a review). With extremely high energies and long traveling distances (> lOMpc), 
the flavor ratio of UHECNs can also help constrain neutrino oscillation parameters [291431] 
and probe exotic physics in the parameter regime inaccessible on Earth. 

To detect UHECNs, enormous amount of matter is required for the target due to their low 
flux and tiny interaction cross section. There are four major detection strategies depending 
on either neutrinos interact with nucleons via the neutral current (NC) interaction, vi + N ^ 
vi + X, or via the charged current (CC) interaction, ui + N — t- + X, where I stands 
for lepton and X for hadronic debris that will develop into hadronic showers. The first 
approach is to observe the optical Cherenkov lights emitted by secondary charged particles 
and showers by an array of optical sensors (e.g. photomultiplier tubes) deployed deep in the 
medium, e.g. under-ice arrays such as AMANDA [35J and IceCube |36] at the South Pole, 
and Baikal, ANTARES, NESTOR, NEMO, KM3NET [37] underwater. The second one is 
to detect horizontal or Earth-skimming neutrino-induced air showers, such as Pierre Auger 
Observatory [38] and HiRes [39]. The third approach is to detect acoustic waves generated 
by the showers, which is still in the R&D stage [IQ]. The last and a very promising one is 
to observe the radio Cherenkov emission from the neutrino-induced showers in dense media 
through the Askaryan effect [H] . Showers propagating in dense media would develop about 
20 % of excess negative charges and would emit Cherenkov radiation, which is coherent in 
the radio frequencies up to a few GHz due to the compact shower size. This effect has 
been verified in a series of beam experiments [12]. The radiated power in the coherent 
regime is proportional to the square of net charges, which is roughly proportional to the 
shower energy, making this technique especially sensitive to UHE showers and thus UHECNs. 
Another advantage of this approach is the long radio attenuation length in some natural 
media, e.g. the Polar ice and salt, with lengths typically of order of 0.1 to 1km, and 
therefore detectors are able to monitor large target volume and achieve greater sensitivity. 
Observatories of this type are: the FORTE satellite [43J looking for neutrino signals from 
the Greenland ice; the balloon-borne antenna array ANITA [:44J overlooking the Antarctic 
ice; and radio telescopes looking for signals from the lunar regolith, e.g. GLUE [45J and 
^ In fact, most of these references consider neutrinos with energy > PeV. 
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LUNASKA [16]. There are also attempts to deploy antennas inside the target media in 
order to lower the threshold energy (to about 10^^ eV), e.g. under-ice antenna arrays RICE 
HH, ARA [M], and ARIANNA gS] in the Antarctic ice; and SalSA [50j in salt dome. 

It is impossible to distinguish between neutrino flavors from the NC interactions be- 
cause their only products are hadronic showers. On the other hand, the charged leptons 
produced in the CC interactions have different energy-loss characteristics in the medium, 
which provides an opportunity to identify the flavor. Optical Cherenkov neutrino telescope, 
e.g. IceCube is able to identify the flavors by event topologies. The muon from the z/^ 
CC interaction leaves a track; whereas showers from all flavor's NC events and z/g's CC events 
lead to localized trigger patterns; and there are double-bang and lollipop events unique to 
Pr induced by the r decays. Beacom et al. [51] proposed a method to deduce the neutrino 
flavor ratio from the measured events of different types, and its feasibility has been widely 
investigated (e.g. [211 Ell EOl E21 E2] ) • However, the instrumented volume, currently of cubic 
kilometer scale for the largest, limits the event rate for UHE neutrinos, which renders it 
challenging to distinguish between z/^ and z/^ events at UHEs as the decay length of r lepton 
exceeds the detector size [53] . 

For radio Cherenkov telescopes such as ARA, the situation is a somewhat different. 
Though this approach cannot detect the track by a single charged particle, the z/^ and 
Pr CC events can in principle be separated through the amount of energy deposited by 
leptons into electromagnetic and hadronic showers [53]. It has also been pointed out that 
different types of showers at UHEs can be distinguished according to their elongation by the 
Landau-Pomeranchuk-Migdal (LPM) effect [511 155]. In addition, z/e CC events can generate 
mixed showers of both types and should have its own characteristic signal feature [HSj- The 
feasibility of this method has been investigated in SalSA [27]. 

Apart from the event signatures, the direction distribution of neutrino events also man- 
ifest themselves in the the energy-loss properties of leptons since neutrinos from different 
directions propagate through different amounts of matter in the Earth. For example, Vt 
can undergo the regeneration process (z/,- — )■ r — )■ u^) without being absorbed due to the r 
decay [531 EH] , and as a result can exhibit a higher flux in the up-going directions. Therefore 
the neutrino distribution can be a useful tool for measuring the flavor ratio and should be 
applicable to any detector with sufficient angular resolution of event direction. A similar 
idea that takes advantage of the event direction distribution has been proposed to constrain 
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neutrino- nucleoli cross sections [59]. 

In this paper we focus on the cosmogenic neutrinos and consider three expected flavor 
ratios when they arrive at the Earth's surface, that is, 1 : 1 : 1 expected in the standard 
oscillation scenario [18j; and 6:1:1 as well as 0:1:1 ratios predicted in the neutrino decay 
models [19] with normal and inverted neutrino mass hierarchy, respectively. The sum of 
ratios is normalized to unity in the following sections. With the detector configuration based 
on ARA [IH] currently under construction at the South Pole, we demonstrate the feasibility 
of extracting the flavor ratio from the event direction distribution, while the inference on 
the flavor ratio at the source is beyond our scope. 

In the next section we present the simulation setup, and derive the expected event direc- 



tion distribution for each flavor in Section III The procedures for hypothetical experiments 
and the extraction of flavor ratios from the direction distribution of pseudo-data are de- 



scribed in Section |IV A| and |IVB[ The successful probability of this method and the flavor 
ratio resolution as functions of the number of observed events and angular resolution of 



neutrino direction are reported in Section IV C 



II. SIMULATION SETUP 



Our simulation code is built by integrating existing packages that consists of two parts. 
The first part is the propagation of the neutrinos and the secondary charged leptons from the 
Earth's surface, where the primary UHECN is generated, to the vicinity of the detector. The 
second part is the simulation of event detection, including the conversion of particle energy 
losses to showers, the development of Cherenkov radiations from showers, the propagation of 
Cherenkov radiations to detectors, and the calculation of detector responses. The planned 
configuration of ARA |18] and the ice properties at the South Pole are adopted as the setting 
in our simulations. 



A. Neutrino Generation and Propagation Using MMC 

The Muon Monte Carlo (MMC) package [60] is employed to generate neutrinos and propa- 
gate all types of neutrinos and secondary charged leptons. In MMC, interaction cross sections 
of neutrinos are evaluated based on Ref . [61] with CTEQ6 parton distribution functions [62] . 
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For charged leptons, energy losses via ionization, pair production, bremsstrahlung, photonu- 
clear interaction, and decay are taken into account, so is z/,- — t- r — )■ v^. regeneration. The 
Kelner-Kokouhn-Petrukhin (KKP) [53] and Bezrukov-Bugaev (BB) [01] parameterizations 
are chosen for cross section calculations of bremsstrahlung and photonuclear interaction, 
respectively. We do not propagate secondary electrons and regard them as losing all of their 
energies to shower developments within a short distance once they are generated. 

Monoenergetic neutrinos are generated isotropically at the Earth's surface and start their 
propagation to the detection volume. The Earth model provided in MMC simulation code 
is used, where the density is calculated based on the Preliminary Reference Earth Model 
(PREM) [65] while the composition as well as the topography on the Earth's crust are not 
considered. The detection volume, which is the South Polar ice sheet in the vicinity of ARA, 
is approximated by a cylindrical ice volume, centered at 1 km below the ice surface with a 
radius of 8 km and a height of 2 km. In the propagation, neutrinos and secondary charged 
leptons are traced until they are either absorbed by the Earth or exiting the detection 
volume; energy losses greater than 1 PeV are treated stochastically. Only those neutrino 
events traversing the detection volume with at least one interaction and with energy loss 
exceeding 1 PeV are reserved for the event detection. The cutoff energy is chosen based on 
the consideration that the radio Cherenkov signals emitted by showers below this energy 
would not be strong enough to trigger the detector efficiently. 



B. Neutrino Event Detection Using SADE 

The Simulation of Askaryan Detection and Events (SADE) package [66] is used for sim- 
ulating neutrino detection. Neutrino events recorded in the previous step, as described 



in Section |II A| are processed individually. Every energy loss exceeding 1 PeV within the 
detection volume is converted into showers of corresponding types. Hadronic products of 
neutrino CC and NC interactions, r decay, and photonuclear interactions are turned into 
hadronic showers, while secondary electron, pair production, and bremsstrahlung are turned 
into electromagnetic showers. In general, a neutrino event can generate multiple showers 
in the detection volume. For example, a Ve event having CC interaction would generate a 
hadronic shower and an electromagnetic shower, whereas a z/^ or a Vr event may have several 
to even hundreds of showers produced by the secondary /x or r lepton, respectively. 



In SADE, shower characteristics and the frequency spectrum of Cherenkov radiation 
are calculated with analytic formulae, where the latter is primarily based on Ref. [07] but 
with a little bit modification of parametrization to account for the decoherence due to the 
longitudinal and the lateral spreads of shower (see Appendix |A] for details). 

A ray-tracing routine then finds the paths of both the direct and the reflected (due to 
the ice-air interface at the surface) rays connecting each shower to each antenna in the ice 
whose index of refraction varies with depth. The flight time, the radiation spectrum taken 
into account the frequency-dependent attenuation along the path as well as the polarization 
at the antenna are calculated for each ray. The spectrum is then Fourier-transformed to the 
electric field received by the antenna in the time domain accordingly. 

Ice properties such as the index of refraction [68], the temperature |69j, and the radio 
attenuation length ^81 170] are based on the results of in situ measurements at the South 
Pole. The index of refraction n and the ice temperature T (in °C) depend only on the depth 
(in km), \z\, 

n{z) = 1.78 - (1.78 - 1.35) exp(-13.2|2|), (1) 
T{z) = -51.5 - 0.45319|z| + 5.822|^|2. (2) 

The attenuation length in turn depends on the ice temperature as well as the radiation 
frequency, and is plotted in Fig. [1} Note that the rising ice temperature with increasing 
depth renders the attenuation length shorter and thus suppresses the detectability of signals 
originating from the bottom part of ice. The birefringence of South Polar ice, which is the 
polarization dependence of the wave speed and the attenuation due to the crystal anisotropy 
and orientation of ice, is not considered here. It is reported [71] that the birefringence is 
observed at the bottom half of the ice sheet and will reduce about 5 % of the neutrino 
detection volume. 

The detector based on the planned configuration of ARA [18| is a hexagonal array of 
37 antenna stations arranged in a triangular grid with 2 km spacing. The array covers a 
total area of about 100 km^. Each station is an autonomously operating cluster of eight 
vertically polarized (Vpol) and eight horizontally polarized (Hpol) antennas evenly deployed 
on four vertical strings with each string placed on one vertex of a square. Each string is at 
a maximum depth of 200 m and is loaded with two antenna pairs, where each antenna pair 
contains a Vpol antenna and an Hpol antenna. The values of parameters for the detector 
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FIG. 1. Attenuation length used in SADE [? ] as a function of the radiation frequency at different 
temperatures (— 30°C, — 40°C, — 45°C and — 50°C, top panel) and as a function of the depth at 
different frequencies (300 MHz, 500MHz and 700MHz, bottom panel), respectively. 
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settings are listed in Table |Tj 

Given the incident electric field E at an antenna, the received signal voltage of the antenna 

Signal is 

Kignal ~ 2^ ' 

with 



where hes and are the effective height and the effective area of the antenna, respectively; 
Zant the antenna impedance, Zq ~ 377Q the impedance of free space, n the index of refraction 
of surrounding medium, G the antenna gain, and c the speed of light in vacuum. The 
direction of effective height for Vpol antennas is in the vertical direction z, while it is in the 
azimuthal direction for the Hpol. In the simulation, the frequency response of the antenna 
is assumed to be a single perfect passband, and the effective height is approximated by a 
single value evaluated at the central frequency of the passband. For each antenna, signals 
coming from different showers are summed in the time domain. We neglect the contribution 
from the reflected signals because they would have sufficient time delays and would suffer 
more attenuation than the direct ones due to longer path length and are more difficult to 
be reconstructed after passing through the less compact snow near the surface (firn) where 
the index of refraction changes rapidly (see Eqn. [T]). 

The root mean square (RMS) thermal noise voltage of an antenna Kms is defined as 

Vrms = ^/^B^sys-^ant-S, (6) 

where is Boltzmann's constant, T^yg the system noise temperature of the receiving antenna 
system, / the radiation frequency and B the frequency bandwidth of the antenna. The values 
of antenna parameters used in the simulation are summarized in Table |T} 

The trigger conditions for a detected event require that i) the received voltage of a 
triggered antenna should exceed three times of its RMS noise voltage (i.e., Vsignai > 3yrms); 
a) at least five out of sixteen antennas in a station are triggered; and 
Hi) at least one station is triggered. 



Monoenergetic neutrinos with initial energies logiQ^Ey/eV) = 17, 17.5, 18, 18.5, 19 and 
19.5 are generated separately, and the numbers of events generated and detected are listed 
in Table [ll] (see Appendix) . 

TABLE I. Parameters for the detector configuration and the antennas, and their values used in 



the simulation. 

Parameter (unit) Value 

Station spacing (km) 2 

Radius of string (m) 10 

Number of strings per station 4 

Separation between paired antennas (m) 5 

Vertical spacing between antenna pairs (m) 20 

Maximum antenna depth (m) 200 

Number of Vpol antennas per station 8 

Number of Hpol antennas per station 8 

Vertical antenna configuration V,H above V,H 

Vpol antenna bandwidth: By (MHz) 150-850 

Hpol antenna bandwidth: Bu (MHz) 200-850 

Antenna impedance: Zant {^) 50 

Antenna gain: G 1.64 

Effective height of Vpol antenna: |/ieff,v| (cm) 11.8 

Effective height of Hpol antenna: |/iefr,H| (cm) 11.3 

System noise temperature of antenna: Tsys(K) 300 

RMS noise voltage of Vpol: Vrms,Y (V) 1.20 x 10^ 

RMS noise voltage of Hpol: Frms.H (V) 1.16 x 10"^ 

Antenna trigger threshold (14ms) 3 

Station trigger threshold (antennas) 5 



10 



III. ANGULAR DISTRIBUTION OF INCOMING NEUTRINOS 



To acquire the expected direction distribution of detected neutrino events with initial 
energy and flavor a, Da{E^, cos 6), the following information are needed: 

a) the flux of isotropic cosmogenic neutrinos for all flavors, ^i,{E„), as well as the incident 
flux ratio among three flavors at the Earth's surface, fe'-fj^'- fr i with = 1; 

b) the interaction probability Pint,a{Eu,cos6), which is defined as the probability that 
a neutrino or its secondary lepton traverses the Earth in the zenith direction cos 6' with- 
out being stopped and interacts (with at least one energy loss exceeding 1 PeV) inside the 
detection volume, to account for the propagation effect; and 

c) the detection efficiency to the subsequent shower(s) generated by the neutrino event 
interacting inside the detection volume, edet,a{,E^,cos9). That is, 

D4E,,cose) = y^^,{E,)PintAEu,cose) (7) 

Xedet,a(^i/,C0S6'), 

with the normalization factor 

N = J2jdEuj dcos9f^^,{E,) (8) 
xPintAEu,cos6)edetAEu,cos9), 

where the distribution has been normalized as a probability distribution independent of the 
total event rate, a = e, fi, r, the superscript E indicates quantities on the Earth, and 9 is 
the angle between the local vertical axis of the detector (z) and the direction of neutrino 
momentum. Hence, events with negative cos 6 are down-going, whereas those with positive 
values are up-going. The neutrino energy range considered in this article is logiQ^E^/eV) = 
16.75-19.75. 



A. Flux and Flavor Ratio on Earth 

The cosmogenic neutrino fluxes for all flavors, ^,y{E,y), have been theoretically predicted 
in Refs. p3l 1721 173] (see Fig. [2]). In the following analysis, the neutrino flux from Ref. [72] 
(red solid curve; hereafter, ESS) is assumed. Note that it is the spectral shape that affects 
the event direction distribution instead of the overall flux normalization. Therefore one 
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FIG. 2. Differential cosmogenic neutrino fluxes predicted by [75] (ESS, red solid curve), [7H] (green 
dashed), the optimistic (purple dashed) and the plausible (blue dashed) scenarios in [H]. The 
shaded region indicates the neutrino energy range considered in the analysis, logj^o(-^i'/6V) = 16.75- 
19.75. 

expects that neutrino fluxes predicted in [72], [73] (green dashed curve), and the optimistic 
scenario in [13] (purple dashed curve) should yield similar distributions. The flux predicted 
in the plausible scenario in [T3] (blue dashed curve) has a steeper spectrum than others, and 



we will present its results later in Sec. IV 



The flavor ratio of neutrinos arriving at the Earth's surface adopted in our analysis is 
fe '■ '■ fr = 1/3 '■ 1/3 : 1/3, as expected in the standard oscillation scenario [T8] : 
0.75 : 0.125 : 0.125 and : 0.5 : 0.5 predicted in the neutrino decay scenarios with normal 
and inverted neutrino mass hierarchy, respectively [19]. Throughout this paper these flavor 
ratios are assumed to be energy-independent over the considered neutrino energy range. 
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B. Interaction Probability 



To obtain the interaction probability Pint from the simulation results, we first divide the 
zenith angle of neutrinos cos 6* into 100 bins with a width of 0.02. The probability at each 
bin is defined as the ratio of the number of survival events inside the detection volume 
to the number of initial incoming events at the Earth's surface. The probability does not 
depend on the azimuthal angle of the neutrino due to the axial symmetry of our Earth model 
and detection volume. The probabilities Pint,a(-E',/, cos 6^) for initial neutrino energies Ei, = 
10^^ eV, 10^^ eV and 10^^ eV are shown in Fig. g 

For different flavors of neutrinos with the same initial energy, the interaction probabilities 
are about the same as cos^^ approaches —1 where neutrinos impinging directly downward 
into the detection volume, whose size is much smaller than the neutrino interaction length 
(about several hundred km [M]). Thus the probability is approximately equal to the size of 
the detection volume divided by the neutrino interaction length for energy transfers greater 
than 1 PeV. The probability for neutrinos impinging downward increases with neutrino 
energy as the neutrino interaction length decreases. 

The probabilities then increase with cos 6' and reach a maximum near the horizontal 
direction (cos 6* ~ 0), where the traveling distance of neutrinos becomes comparable to the 
interaction length and the detection volume has its maximum span. The probabilities for 
different flavors diverge due to the difference in the energy-loss property between different 
flavors of charged leptons produced in the CC interactions. The longer the lepton can 
propagate, the higher the probability is. Contrary to electrons, which would be stopped 
immediately after their creation and would develop into electromagnetic showers, EeV n 
and r leptons can on the average propagate distance of order of 10 km before come to a stop 
[75] . Muons lose their energy mostly via pair production, while tauons via pair production 
as well as photonuclear interaction. So the probabilities of flnding z/^ and u^- are higher 
than that for z/g. The interaction probability in these directions also increases with neutrino 
energies because of the decrease of the neutrino interaction length and the lepton propagation 
range. 

The probability for up-going neutrinos (cos > 0) is suppressed as the neutrino traveling 
distance becomes longer than the interaction length. The higher the initial neutrino energy, 
the shorter the interaction length, and hence the distribution terminates at smaller cos^. 
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The Earth attenuates the neutrino both in energy through NC and CC interactions and in 
number through the stoppage of the charged lepton produced in the CC interaction. As a 
special case, r lepton, having a decay length of about 50 x (i?^/PeV) m, can transform to 
Vr through decay before losing too much energy |53l EH]. Therefore Vr coming from below 
the horizon would have an apparent larger probability than other two flavors. This is a 
critical feature for the flavor ratio determination proposed in this article, as we will further 
ellaborate below. 

C. Detection Efficiency 

The detection efficiency edet defined here depends not only on the nature of Cherenkov 
radiation, the ice properties and the detector configuration, but also on how neutrinos and 
secondary leptons deposit their energies in the shower. But for the simplicity of compu- 
tation, we do not decompose it further into a product of the probability that a neutrino 
or a secondary lepton in specific direction generates shower(s) of specific energy at specific 
position, times the detection efficiency to each individual shower. 

Similar to the definition of the interaction probability, the detection efficiency is defined 
as the number of events detected divided by the number of events interacting inside the 
detection volume. The detection efficiency is hence an averaged quantity over the azimuthal 
direction and the shower position. The results for different initial neutrino energies and 
fiavors are shown in Fig. |4j 

The efficiency increases with neutrino energy simply because shower size is proportional 
to neutrino energy. The Cherenkov cone generated by neutrinos events with cos^^ ^ —0.7 
travels downward and has an opening angle 56° in ice, so it is less possible to cover the 
area where the antennas are deployed. This leads to the common cutoff at cos^^ ^ —0.7 
for all neutrino energies. In Fig. |4| we see the fluctuations of the efficiencies for z/e (top 
panel) and z/^ (middle panel) at cos^ > 0.1. This is due to the smallness of the number of 
events arriving at the detection volume. Such results are therefore not reliable. However 
one generic feature remains valid; that is, the up-going neutrino events diminish since their 
energies are severely damped by the Earth. 

The detection efficiency for is the highest among the three fiavors, because once CC 
interaction occurs all the neutrino energy is released into showers and strong signals are 
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FIG. 3. Interaction probability Pint (see text for definition) as a function of the neutrino zenith 
direction cos^ for fg (red dash-dotted), (green solid), and v-r- (blue dashed), and for initial 
energies Ei, = 10^^ eV, lO^'^eV and 10^^ eV (from top to bottom panel). 



emitted. For z/^, although there are plenty of electromagnetic showers produced by muon, 
these showers tend to have lower energies so that they are less likely to be detected by the 
sparse antenna array. This leads to lower detection efficiency of z/^, where the NC-induced 
hadronic showers account for about 80 % of detected events for 1 EeV v^. The situation is 
similar for v^, where most hadronic showers from photonuclear interaction and electromag- 
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netic showers from pair production do not trigger detector efficiently while hadronic showers 
induced by NC interaction and r decay (r — )■ z/,- + hadrons) account for the most detected 
events. The convergence of efficiency for up-going z/^-s of different initial energies results 
from the degradation of neutrino energy to few PeV by the regeneration process and NC 
interaction |53] . 

D. Event Direction Distribution and All-Sky Flavor Ratios of Events 

Because the initial neutrino energy is sampled only with discrete values of equal logarith- 
mic interval A = Alogio-E'i/ = 0.5 in the calculation of Pint and edet, the direction distribu- 
tion integrated over the j-th energy bin ranging from log^^gi^j — A/2 to log^^g-^j + ^/2, is 
approximated by 



with log;^Q(-E'j/eV) = 17, 17.5, 18, 18.5, 19 and 19.5. The results assuming the ESS neutrino 



been normalized to unity and the relative fraction contributed by each bin is also shown. 
It appears that EeV neutrinos contribute the most to the detected events for every flavor, 
because of the compromise between two competing effects: the decrease in the neutrino flux 
versus the increase in the detection efficiency with neutrino energy. 

Finally, after summing over all distributions of different energy bins, the expected event 
direction distribution for each flavor Da{cos6) is obtained, which is shown in Fig. [6] We 
define the "all-sky" flavor ratio of detected events, fe '■ fp. '■ fr, as the event ratios among 
flavors after integrating over all zenith directions. If the ESS neutrino flux and incident 
flavor ratios at the Earth's surface, : : = 1/3 : 1/3 : 1/3, are assumed, we flnd 
that fe '■ ffi '■ fr = 0.584 : 0.153 : 0.263. The z/g events account for the most portion because 
of their higher detection efficiency. The Uj. events have a different shape compared to the 
other two flavors especially in the horizontal and up-going directions, primarily due to its 
special interaction probability (see the graph in the right panel. Fig. |6]). These will be used 
to extract the flavor composition at the Earth's surface (/^s) in the next section. However, 
the resemblance between z/g and z/^ distributions will lead to the degeneracy in the flavor 




X Pint,a {Ej , COS 6^) edet,a {Ej , COS 6*) , 




flux are plotted in Fig. [5} where the total area under the distributions for each flavor has 
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FIG. 4. The detection efficiency edet (see text for definition) as a function of neutrino zenithal 
direction cos^ for initial energies logj^gl-^f/sV) = 17 (filled red squares; red solid line), 17.5 (open 
red; red dashed), 18 (filled green; green solid), 18.5 (open green; green dashed), 19 (filled blue; 
blue solid), 19.5 (open blue; blue dashed), and for flavor u^, v^, and Ur (from top to bottom panel, 
respectively) . 



ratio extraction, and an extra constraint is required, for example, the v^-Vr symmetry. 

Once given reference incident flux ratios /^g ■ //fo • frfl ^^cl the corresponding event 
ratios /e,o : //i,o • /t.Oj the all-sky event ratios fe '■ ffi '■ fr for arbitrary incident flux ratios 
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fe '■ fl^ '■ fr can be deduced as 

f fE fE £E 

~ ~f fE fE fE < f fE fE fE If fE fE fE^ ^ ^ 

JeflJe J fiflJrfl "r J tJ.flJ eflJ J rfl "r JrflJeflJ fiflJr 

and similarly for and fr- The conversion from /s to f^s can be done by just interchanging 
/ with f^. The relation between z/g event ratio fe and flux ratio is plotted in Fig. [7| 
assuming ESS neutrino flux. 

IV. PSEUDO-OBSERVATION AND FLAVOR RATIO EXTRACTION 

To investigate the discriminating power of flavor ratio reconstruction using event direction 
distribution, we generate pseudo-observation data from simulated events and fit the direction 
distribution to extract the flavor ratios, and repeat the processes to determine the statistical 
uncertainty of the extracted ratio. Results assuming different incident flux ratios, numbers 
of observed events, and angular resolution of detector are then presented. 

A. Pseudo-Data Samples 

The pseudo-data sample is prepared in three steps. First, detected neutrino events are 
randomly picked out one by one from the database of simulated events according to the 
expected direction distributions, Da{cos6), until any of direction bins of the database ex- 
hausts its events. Secondly, the zenith angle, 9, of each pick-up event is smeared by adding a 
Gaussian distributed random number with zero mean and the standard deviation A6 equal 
to an assigned experimental angular error in reconstructed neutrino direction. Finally, this 
event collection is evenly divided into subsets, where each data set represents a hypothetical 
experimental data sample with a total number of detected events equal to Nohs and these 
data sets constitute a statistical ensemble. 

In the following, A^^obs varies from 50 to 500 with an increment of 50, while A^^ from 0° 
to 6°. 

B. Fitting Pseudo-Data 

To construct the fitting function for the hypothetical experimental data samples with 
angular resolution A^, the expected direction distribution for each flavor a is convolved 
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FIG. 5. The expected direction distribution integrated over the energy bins, log^oC-^j) ^ 0.25, 
of initial neutrinos, where log;^Q(£'j/eV) = 17 (red solid), 17.5 (red dashed), 18 (green solid), 
18.5 (green dashed), 19 (blue solid), 19.5 (blue dashed), for flavor fg, i^^i, and Ur (from top to 
bottom panel, respectively). The distributions have been normalized so that each one represents 
its fractional contribution to the corresponding flavor, and the relative fraction of each energy bin 
is listed in the legend. The ESS neutrino spectrum [72] is assumed. 
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FIG. 6. Left panel: the expected direction distribution of detected neutrino events for each flavor 
a, Da{cos9), assuming the ESS flux [72] and the incident flux ratios /f^ : : = 1/3 : 1/3 : 1/3. 
The all-sky flavor ratio of detected events fe ■ ffi '■ fr is 0.584 : 0.153 : 0.263. Right panel: the 
expected direction distributions are scaled to equal fraction, 1/3 for each flavor, for comparing the 
curve shapes. 
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FIG. 7. The relation of the fraction of detected Vf, events, /e versus the Vf. fraction of incident 
neutrino flux at the Earth's surface, /^f, where the ESS flux model f72] is assumed. The values for 
the standard case (/^f :/^:/^ = 1/3:1/3:1/3, blue squares) and for the neutrino decay case 
(0.75 : 0.125 : 0.125, red stars) predicted in [19] are labeled on the axes. 

with a Gaussian resolution function of standard deviation in Q space, G{6';6, (Ad)'^). 
The convolved distribution, 

Ma{e;Ae)(x j D^{e)G{e';e,{Aef)de\ (ii) 

is then normalized so that 




Ma{cos9]Ae)dcos9 = 1, 
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and the probability density function of event direction can be expressed as 



P{cos9; /i, /2, A9) = /iMe(cos0) + (1 - A) (12) 

x[/2M^(cos^) + (1 - /2)M.(cos^)], 

where /i and /2 are unknown fraction coefficients with values between zero and one. This 
expression ensures /i and /2 are independent of each other, and they are simply related to 
the event flavor ratios by fe = fi, = (1 - /i)/2, and = (1 - - /2). 

The maximum likelihood estimation is applied for flavor ratio extraction. For an experi- 
mental data set with total events of A^^obs, the likelihood function is deflned as 

L(cos^,;/i,/2,A^) = n^(cos^.;/i,/2,A^), (13) 

i=l 

where the subscript i = 1,2, . . . , A^obs stands for the i-th event. The true values of /i and /2 
are estimated by maximizing L, or equivalently minimize the negative log-likelihood (NLL), 

- In L = - ^ In P(cos 9,; f,, A9) . (14) 

2=1 

We perform grid search in the parameter space to flnd the minimum of NLL and the associ- 
ated best flt values of /i and /2, denoted as fi (or /g) and /2. A reconstruction is deflned as 
failed if /i or /2 is outside either boundary (zero or one). It may result from the insufficient 
statistics of events, or that the true value is near the boundary. 

However, an extra constraint is required to avoid multiple solutions arising from the 



similarity in distribution shape between z/g and z/„, as pointed out in Section HID The z/^- 



z/t- symmetry, i.e. the incident z/^ and u^. fluxes on the Earth are of equal amount (/^ = /^), 
is imposed for this purpose. As a consequence, the value of /2 is known and flxed, and there 
is only one parameter f\, the fraction of u^. events, left to be fltted. In addition, for simplicity, 
we assume the observer knows the exact shape of neutrino spectrum; that is, the neutrino 
fluxes assumed for generating pseudo-data samples and fltting function are identical. The 
latter constraint will be relaxed in the next section and the mismatch in shape between real 
and expected spectra will introduce a systematic bias to the extracted ratio. 



When doing the fitting, /i and /2 arc permitted to have unphysical values. 
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C. Successful Probability and Flavor Ratio Resolution 

After every data set in the ensemble has been fitted for a given A'^obs and A^, the successful 
probability of flavor ratio extraction is calculated as the number of successful data sets 



indicate the uncertainties due to finite number of data sets and the shortest 68.3 % Bayesian 
credible intervals with uniform prior are used (see Appendix [b| or [76] for details). |^ The 
probability for initial ratio of : 0.5 : 0.5 is not shown here because it is always around 50 % 
regardless of the number of observed events and angular resolution assumed. 

In general, the successful probability increases with the number of observed events and a 
better angular resolution, as expected. The probability for the standard scenario (left panel) 
is greater than those for the decay scenarios (the one with normal hierarchy is shown in the 
right panel) since the latter have true z/g ratios so closed to the boundaries that even small 
statistical fluctuation may result in failure. 

For the standard scenario, the successful probability is greater than 70 % for Aobs > 100 
and becomes over 90 % for Aobs > 200, if the angular resolution is within 6°. On the 
other hand, for the decay scenario with normal hierarchy, the probability is over 50 % once 
Aobs > 100. It was reported in Ref. |l8] that fully deployed ARA is able to detect about 50 
cosmogenic neutrinos in three years. Therefore our method is feasible for ARA to extract 
the flavor ratio of cosmogenic neutrinos. 

We define the resolution of event fraction as the spread of all fitted values /eS in the 
ensemble with respect to the expected value /e,exp, 



where Ag is the number of hypothetical experimental data sets. If the fitting is not disturbed 



distribution as the number of observed events Aobs increases, and therefore the resolution 
defined here is approximately Icr uncertainty at 68 % confidence level. In order to reduce 
the effect induced by boundary cutoff, the resolution is calculated separately at both sides 
of the expected value and denoted as -R+ and -R_. 

^ The concept of credible level is somewhat different from that of confidence level, because for the former 
our knowledge (the prior distribution) has been imposed. 



divided by the total number of data sets. The result is shown in Fig. [8} where the error bars 




(15) 



by the boundary cutoff, then the distribution of /eS asymptotically approaches a Gaussian 
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FIG. 8. The successful probability of flavor ratio extraction as a function of the number of 
observed events A^^^bs, assuming initial ratio of 1/3 : 1/3 : 1/3 (left panel) and 0.75 : 0.125 : 0.125 
(right panel), where A'^obs ranges from 50 to 500 with an interval of 50. Results for different 
angular resolutions are plotted: = 0° (red), 2° (green), 4° (blue), 6° (black). The error bars, 
which result from the statistical uncertainty due to the finite number of data sets in the ensemble, 
represent the shortest 68.3% Bayesian credible intervals with uniform prior [76] (see Appendix [B] 
for details). Each curve is horizontally shifte ghtly for clarity. The ESS neutrino flux [72] is 
assumed in both data generating and fitting, and the v^-Vt symmetry is assumed in fitting so that 
the relative ratio between them is known and fixed. 



The upper and lower bounds of z/e event fraction, /e,exp ± R±, are then transformed into 
the corresponding flavor ratios at the Earth's surface, f^tme ^ -^±' according to Eqn. 
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(with / interchanging with f^; see also Fig. [t]). The result for given A^^obs and is shown 
in Fig. [9} The upper bounds for the initial ratio of 0.75 : 0.125 : 0.125 (blue error bars) are 
less meaningful because of the boundary cutoff. The resolution of Ue ratio on the Earth 



for given Nohs and is plotted in Fig. 10 The lower (upper) resolution is taken in the 
decay scenario with normal (inverted) hierarchy. The resolution in the standard scenario 
is asymmetric on both sides, with the upper resolution worse than the lower one typically 
by about 0.06, due to the nonlinear conversion relation between event ratios and flux ratios 
(Fig. [7]), and hence the average value is taken. 

To measure the separability of a flavor ratio predicted in a "fake" scenario from that 
extracted from a given "true" scenario, we define the discriminating power of flavor ratios 
as 

J. . • /e, model /e,exp /1 /^\ 

discrimmatmg power = [lb) 

R 

where /e,modei is the event ratio expected by the scenario to be examined, and /e,exp and R 
are the expected u^. event ratio and its resolution in the assumed true scenario, respectively. 
Note it is the z/g "event" ratio and its resolution that are used in the definition because of 
the asymptotic normality of its distribution, whereas the distribution of Ug flux ratio is 
skew due to the nonlinear conversion relation. The discriminating powers assuming different 



true scenarios are shown in Fig. 11 



In Ref . [H] , the ARA simulation result reports an angular resolution of neutrino direction 
about 6°. Under this circumstances, given 100 observed neutrino events, the successful 
probability is about 75 % in the standard scenario while about 50 % for both decay scenarios 
(black curves in Fig. |8]for the one with normal hierarchy). If the standard scenario is the 
real one, the discriminating power with respect to the decay scenario with normal (inverted) 
hierarchy will be at about l.lcr {2.2a) level. If the decay scenario with normal (inverted) 
hierarchy is real, then the discriminating power with respect to the standard scenario will be 
at about 0.8cr (1.4(t) level, while that to its inverted (normal) counterpart will be at about 
2.4(7 (2.2(t) level. Hence a preliminary constraint on the neutrino decay can be set by the 
method proposed in this article. 

If the number of the observed events is doubled with the angular resolution flxed (A'obs = 
200, A9 = 6°), the successful probability rises to 90% for the standard scenario while still 

25 



around 50 % for the decay ones; the discriminating power between the standard scenario and 
either of the decay scenario increases by about 0.2a while by about OAa between the decay 
scenarios. If the angular resolution is improved by 4° with the number of observed events 
fixed {Nohs = 100, A9 = 2°), the successful probability becomes about 83 % for the standard 
scenario while about 55 % for the decay one with normal hierarchy; and the discriminating 
power between the standard scenario and the decay one with normal hierarchy increases by 
about 0.2cr, and by about O.lcr for the decay one with inverted hierarchy to the other two 
while by about OAa in the contrary case. The limit of our method with perfect angular 



resolution is also shown in Fig. [8| Fig. 10 and Fig. 11 (red lines). 

To accumulate more neutrino events, apart from waiting for more neutrinos to come, one 
can increase the event rate by extending the antenna array. To achieve a higher angular 
resolution of neutrino direction, on the other hand, an antenna array with denser grid is 
required for a better imaging of the Cherenkov cone. Therefore the results presented here 
can serve as a reference for optimizing the future detector configuration of ARA or other 
similar observatories. 



D. Systematic Uncertainty from Neutrino Spectrum 



Fig. 12 , Fig. 13 , Fig. 14 and Fig. [15] show the results assuming the neutrino flux predicted 
in the "plausible" scenario in Ref. [H] (blue dashed curve in Fig. [2]). The successful proba- 
bility is lower than that assuming the ESS flux and the resolution is worse. Recall that the 
flux in Ref. [14j is steeper than the ESS one, i.e. there are more neutrino events with lower 
energies. Since the interaction probabilities between flavors are less distinct from each other 
at lower energies (see Fig. [3]), the event direction distributions between flavors are harder to 
be distinguished. 

If the observer does not exactly know the "real" neutrino spectrum (i.e., the one used 
to generate pseudo-data), then the mismatch in shape between the real and the expected 
(i.e., the one used to fit the data) spectra will introduce a systematic bias to the extracted 
ratio. We found that if the real spectrum is the plausible scenario in [14] while the spectrum 
expected by the observer is the ESS model, the bias is about 10 %. More specifically, the 
average fitted value drops from 0.58 to 0.53, because a steeper neutrino spectrum leads to 
more events coming from below the horizon due to the longer interaction length of neutrinos 
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FIG. 9. The resolution of Ve ratio at the Earth's surface plus the true ratio ff^R^ (see the text for 
definition) versus the number of observed events A^'obs for initial ratios /^f:/^:/.^ = 1/3:1/3:1/3 
(green, middle), 0.75 : 0.125 : 0.125 (blue, top), and : 0.5 : 0.5 (red, bottom), where A'^obs ranges 
from 50 to 500 with interval of 50. Results assuming different angular resolutions are shown: = 
0° (the perfect case, top left panel), 2° (top right), 4° (bottom left), and 6° (bottom right). The 
error bars for different initial ratios are horizontally shifted slightly for clarity. The ESS neutrino 
flux [72] is assumed in both data generating and fitting, and the v^-Vr symmetry is assumed in 
fitting so that the relative ratio between them is known and fixed. 



with lower energies. In order to fit the event distribution from the steeper spectrum with 
the flatter one, the Vr fraction has to be increased while the z/g fraction decreased. 



E. Systematic Uncertainty from Neutrino Cross Sections 

The calculation of neutrino cross section in Ref. [61] reports an uncertainty of factor of 
2^^ at around 10^° eV. Therefore we vary the neutrino cross section by a factor of two, 

27 



Ratio: 1:1:1, Station spacing: 2km 



n 
ui 

c 0.4 
□ 

o 

re 



01 

> 



0.35 



5 0.3 



w 0.25 
DC 



0.2 
0.15 



l_ I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I M 

Z_ ESS2001 Jl 

~ m ~ 

-&A9 = 0.0° ♦A 9 = 2.0° 

-s-A9 = 4.0° ■•■A 9 = 6.0° 




I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 



I I I I I I I I I I I - 



50 



100 150 200 250 300 350 400 450 500 



Ratio; 6:1:1, Station spacing: 2km 

j 0.65 F 
j 

I I 
: 0.6- 
) 

I 0.55 E 

:' 0.5 = 
) 

i 0.45 E 

J I 
I 0.4 = 

0.35 E 

0.3 E 

0.25 = 



50 100 150 200 250 300 350 400 450 500 
Number of observed events N 





1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 






ESS2001 






^a6 = o.o° -•-ae = 2.o° 






^AG = 4.0° -•■4Q = 6.0° 




= 












-^.^ \ 
















; 


1 


3 



Number of observed events N 

Ratio: 0:1:1, Station spacing: 2km 

S 0.351" I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 



ESS2001 

^*A9 = 0.0°-»-Ae = 2.0° 
*A0 = 4.0°-. A9 = 6.0° 




50 100 150 200 250 300 350 400 450 500 
Number of observed events N 



FIG. 10. The resolution of Uf. ratio at the Earth's surface (see text for definition) versus the 
number of observed events N^^g for angular resolution of neutrino direction = 0° (the perfect 
case, red solid line and open square), 2° (green dashed line and full square), 4° (blue dash-dotted 
line and open circle), and 6° (black dotted line and full circle). Results assuming initial flavor 
ratios of 1/3 : 1/3 : 1/3 (top panel), 0.75 : 0.125 : 0.125 (middle panel), and : 0.5 : 0.5 (bottom 
panel) are shown. The ESS neutrino flux [72] is assumed in both data generating and fitting, and 
the f^-i^r symmetry is assumed in fitting so that the relative ratio between them is known and 
fixed. 
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Discriminating Power for 1:1:1 to 0:1:1 



Discriminating Power for 1:1:1 to 6:1:1 




FIG. 11. The discriminating power of v^, flavor ratio versus the number of observed events iVobs for 
angular resolution of neutrino direction = 0° (the perfect case, red solid line and open square), 
2° (green dashed line and full square), 4° (blue dash-dotted line and open circle), and 6° (black 
dotted line and full circle). Results assuming initial flavor ratios of 1/3 : 1/3 : 1/3 (top panels), 
0.75 : 0.125 : 0.125 (middle panels), and : 0.5 : 0.5 (bottom panels) are shown. The ESS neutrino 
flux [72] is assumed in both data generating and fitting, and the v^-Vr symmetry is assumed in 
fitting so that the relative ratio between them is known and fixed. 
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FIG. 12. The successful probability of flavor ratio extraction. Similar to Fig. |8] except that the 
assumed neutrino flux is the plausible scenario in Ref. [H] (blue dashed curve in Fig. [2]) . 
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FIG. 13. The resolution of fg ratio on the Earth plus the true ratio. Similar to Fig. [9] except that 
the assumed neutrino flux is the plausible scenario in Ref. [H] (blue dashed curve in Fig. [2]). 

and the results are shown in Fig. 16, Fig. 17 Fig. 18 and Fig. [l9j A larger neutrino 
cross section leads to lower successful probability and resolution of flavor ratio, because the 
shorter interaction length increase the fraction of events in the down-going directions where 
the interaction probabilities for different flavors are more or less the same, while makes the 
more distinct, up-going parts slightly decrease. 



V. CONCLUSION AND FUTURE WORKS 



In summary, measuring the flavor ratio of UHE cosmic neutrinos can not only reveal the 
physical conditions at UHECR sources but also probe neutrino oscillation parameters and 
non-standard physics that might involve during the propagation. Neutrino observatories 
using the radio Cherenkov technique such as ARA [IS], ARIANNA [IH], and SalSA [50], are 

For the results with neutrino cross section divided by two, the number of events in the ensemble are less 
than 5000, so I choose not to show them. 
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FIG. 14. The resolution of Up ratio at the Earth's surface versus the number of observed events. 



Similar to Fig. 10 except that the assumed neutrino flux is the plausible scenario in Ref. [H] (blue 
dashed curve in Fig. [2|). 



sensitive in the UHE regime and expected to accumulate of order of 10 to 100 cosmogenic 
neutrinos per year in the near future, and hence provide sufficient statistics for the flavor 
ratio identification. 

In this work, the direction distribution of neutrino events is proposed to determine the 
flavor ratio. In order to investigate its feasibility, a simulation is constructed and the ex- 
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Discriminating Power for 1:1:1 to 6:1:1 



Discriminating Power for 6:1:1 to 1:1:1 




FIG. 15. The discriminating power of v^, flavor ratio versus the number of observed events. Similar 



to Fig. 11 except that the assumed neutrino flux is the plausible scenario in Ref. [14J (blue dashed 
curve in Fig. [2|) . 



pected distribution for ARA is derived. It is found that Vr distribution has a different shape 
from the other two, but the distributions for and resemble each other. Therefore an 
additional constraint, the z/^-z^r symmetry, is imposed on the fitting of data distribution to 
avoid multiple solutions. 

This method is proved to be feasible for ARA, e.g., for 100 events and an angular reso- 
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FIG. 16. The probability of successful flavor ratio reconstruction. Same as Fig. [8] except that the 
neutrino cross section is multiplied by a factor of two. 
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ratio resolution fStatiort spacing:2km, A 6 = 0.0°) ratio resoiution (Station spacing:2i(m, A 6 = 2.0°) 




FIG. 17. The resolution of v^. ratio on the Earth plus the true ratio. Same as Fig. [9] except that 
the neutrino cross section is multiphed by a factor of two. 

lution of 6°, the successful probability is about 70% in the standard scenario and over 50% 
for the neutrino decay models considered here. This method is also able to set a preliminary 
constraint on the neutrino decay, e.g., given an angular resolution of 6°, it requires about 
250 events to separate the standard scenario from the decay model with inverted hierarchy 
by about la level. Therefore, the flavor resolution as a function of the number of observed 
events and the angular resolution of neutrino direction presented here can serve as a ref- 
erence for optimizing the future conflguration of ARA. Similar procedures can be done for 
other observatories. 

However, we have not fully taken advantage of all the information possessed by neutrino 
events. For example, the spatial distribution and the energies of showers, which can be 
retrieved by incorporating the vertex reconstruction, can help to distinguish electron flavor 
from the other two; and the characteristics of radio signals may allow the classification of 
shower types and further the identification of the neutrino flavor event by event. In addition, 
the near-fleld effect of Cherenkov radiation has to be considered, since in some cases the 
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FIG. 18. The resolution of ratio on the Earth. Same as Fig. 10 except that the neutrino cross 
section is multiphed by a factor of two. 



observation distances are comparable to the shower lengths and the far-field approximation 
applied here would underestimates both the signal strength and the angular width of the 
Cherenkov cone, as pointed out in Ref. [77]. The neutrino spectrum and the cross sections 
should also be incorporated into the fitting in order to reduce the systematic uncertainties. 
These will be included in the future simulation, and we expect to get a higher resolution of 
the neutrino flavor ratio. 
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Discriminating Power for 1:1:1 to 6:1:1 



Discriminating Power for 6:1:1 to 1:1:1 
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FIG. 19. The discriminating power of v^. flavor ratio. Same as Fig. 11 except that the neutrino 
cross section is multiphed by a factor of two. 
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Appendix A: Calculation of Radio Cherenkov Spectrum 

In SADE [? ] the calculation of Cherenkov radiation invokes the far-field approximation 
where the observation distance R is much longer than the longitudinal spread of the shower 
4h, Ish/R ^ 1, and the radiation is assumed to be emitted from the point at the shower 
maximum. The frequency spectrum of Cherenkov radiation at an observation distance R 
and an observation angle 9 with respect to the shower axis is parameterized as 



where v is the radiation frequency, 6 = 6 — 6^ the angular offset from the Cherenkov angle 
n the index of refraction in the medium, /o the normalization factor, sin^ the projection 
factor of the electric field to the line of sight, and y and z denote the transverse direction and 
the longitudinal direction along the shower axis, respectively. The two integrals containing 
the charge distributions of the shower fy and fz (with J dyfy{y) = f dzfz{z) = 1 over 
the shower) times the phase terms with the phases (py and 0^ account for the decoherence 
effect due to the spread of shower along the transverse and the longitudinal directions, 
respectively. If the emission from the entire shower is in-phase, the integrated value is equal 
to one, otherwise it is less than one. The integrals can be further parameterized into 




(Al) 




(A2) 



(A3) 



with the scaled frequencies 



y 



ay{S, n)z/. 



(A4) 



(A5) 



and 



ay{6, n) 




(cos 5 + 



sin 5 tan 6'c 



(A6) 



a^((5,ra) 



(sin 5 + 



tan 6'c tan6'c,o' 
1 — cos 5 , tan Oc 



(AT) 



tan 6c tan 6'c,o 



^ This section is adapted from the technical document of SADE. 
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where Rm is the Mohere radius, the subscript denotes the fiducial value used for scaling 
the spectrum to different medium or longitudinal spread of shower, and the angular terms 
arise from the phase difference. There is a similar discussion in Ref. [78] . 

The value of /o and the function forms of Gy and Gz are based on Ref. |67j . 

z/q sm Bcfl ITeV 

G.K) = exp(-i(^^)2), (AlO) 

with /o in unit of VMHz~^, z/g = 1.15 GHz, 61^0 = arccos(l/1.78), i?sh the shower energy, and 
ciu,AYZ = 4h,o(3.72°/3.1 m). Note that the frequency z/ of the lateral decoherence factor in the 
original parametrization of |67] has been replaced by Uy, and the index of the exponential 
term in the longitudinal decoherence factor has a different expression. 



Appendix B: Uncertainty Estimation of Efficiency 



For the interaction probability and the detection efficiency defined in Sec. Ill, as well as 



the probability of successful fiavor ratio extraction in Sec. |IV[ there is a selection efficiency e 
introduced. Such selection efficiency is defined as the ratio of the number of entries satisfying 
certain condition, k, over the total number of entries, A^: 

k 

with its value between zero and one. For all the three defined efficiencies, /c is a random 
variable due to the underlying stochastic nature of particle interaction process, and whether 
an entry passes or not is independent of each other. Therefore the selection can be regarded 
as a binomial process with given total number of trials and the unknown probability of 
success e, and k follows the binomial distribution. 

Due to the finite number of total entries, the estimated efficiency for given k and is gen- 
erally not the same as the true value, and therefore its deviation should also be disclosed. 
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Since it is shown that neither Poisson error nor binomial error can well represent the uncer- 
tainty when e is near the boundaries [70] , we instead follow the Bayesian approach presented 
in Ref. [TH] to estimate the uncertainty. 

According to Bayes' theorem, the probability distribution of e for given k and N, or the 
posterior probability distribution of e, is 

P{e\k, N) oc P{k\e, N)P{e), (B2) 

where -P(e) is the prior probability distribution of e, which includes our knowledge about 
the parameter e before the data is considered. Because e lies between zero and one, and 
there is no further knowledge about the true value, the uniform prior between zero and one 
is chosen, 

f 1, if < e < 1 
P(e)= - - , (B3) 

I 0, otherwise 

and then the posterior distribution after normalization becomes 

which is a beta distribution with shape parameters of {k + 1) and (A^ — k + 1). The mode, 
i.e. the most likely value, is just k/N a.s expected by intuition, and the shortest credible 
intervals covering 68.3% of the posterior probability are chosen to represent the errors of 
efficiency in this article. 
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